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Abstract 

We present two Monte Carlo algorithms of the Markovian type which solve the 
modified QCD evolution equations at the NLO level. The modifications with respect 
to the standard DGLAP evolution concern the argument of the strong coupling 
constant as- We analyze the z-dependent argument and then the /cr-dependent 
one. The evolution time variable is identified with the rapidity. The two algorithms 
are tested to the 0.05% precision level. We find that the NLO corrections in the 
evolution of parton momentum distributions with /cy-dependent coupling constant 
are of the order of 10 to 20%, and in a small x region even up to 30%, with respect 
to the LO contributions. 
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1 Introduction 



With the first LHC results coming soon the precision of the QCD calculations becomes 
an important issue. As usuall, one finds two approaches to the calculations: fixed-order 
calculations and resumed to all orders Monte Carlo (MC) parton showers. The fixed- 
order results are impressive. Let us give a few examples of such calculations: the splitting 
functions, calculated up to NNLO level [1, 2] or various differential and semi-inclusive 
distributions, calculated in the NNLO and even NNNLO approximation see e.g. [3-5]. On 
the other hand, the MC parton shower approach, indispensable in describing complicated 
experimental signal selection procedures, does not achieve such a high precision. In fact, 
the complete NLO level has not been reached yet by any of the available MC parton 
shower codes. In most of the cases these codes are based on some form of an improved LO 
approximation, see e.g. [6-10]. Another approach is based on combining the NLO matrix 
element with the LO parton shower [11,12]. 

Some time ago some of us started the MC parton shower project with the ambitious 
goal of reaching the NLO precision level. The project is based on the DGLAP evolution 
equation [13] at the NLO level [14-17]. We began by creating two different algorithms, 
and then constructing two MC programs, that solve the evolution equation and simulate 
the one-hemisphere collinear parton shower. First algorithm, called EvolFMC (Markovian 
Monte Carlo), is based on the principle of a Markovian process. It generates the DGLAP- 
type evolution of Parton Distribution Functions (PDFs) up to the NLO level including 
both gluons and quarks. In Ref. [18] we have shown for the first time that with the 
aid of modern CPU such a MC code can solve the LO evolution equations with the 
precision below 10~ 3 , i.e. comparable or even better than the other numerical methods. 
This study has been extended to the NLO evolution in Ref. [19]. The second algorithm 
called CMC is an entirely new type of an algorithm, an example of a wider class that we 
named " Constrained Markovian Monte Carlo" . CMC is designed to solve the Markovian 
evolution with superimposed constraint on both x and flavor type of the outgoing parton 
[20,21]. Such an evolution with predefined end-point, mandatory for the simulation of the 
resonant processes, is an alternative for the commonly used workaround solution of this 
problem based on the so-called "backward evolution" [22]. The ultimate goal is to combine 
two initial-state non-collinear NLO constrained evolutions (in two hemispheres) together 
with the hard process into one NLO parton shower for the Drell-Yan-type processes, either 
in full agreement with the DGLAP evolution or including effects of a modified-DGLAP 
type. 

The EvolFMC code includes also the modified DGLAP evolutions in which the argument 
of the coupling constant is replaced by more complicated functions of x, x' and Q 2 [23-27]. 
In particular, the use of kr as the argument is known to effectively resum some of the 
higher order soft corrections [23,26,28]. For this reason it is the preferred choice of most 
of the MC parton shower codes. The is used as an argument also by the CCFM 
evolution equation [29]. This equation is designed to effectively interpolate between the 
DGLAP evolution in the large x region and the BFKL-type behavior in the low x region. 
The CCFM equation uses, apart from the modified arguments of coupling constant and 
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angular ordering, also the "non-Sudakov form factor", important in the small x region. 

In the case of the developed by us EvolFMC code the option of modified argument of 
the coupling constant has been so far implemented only in the LO approximation [30,31]. 
This is a significant limitation since both the NLO kernels and the /^-dependent coupling 
constant are important for the final precision of the parton shower. In the present paper 
we will fill in this gap and present the NLO extension of the /c^-dependent EvolFMC code. 
In this scheme the evolution variable is understood to be the rapidity of emitted partons, 
ensuring the angular ordering. We will discuss also another choice of the argument: 
Q(l — z). It is based on the variant of the CCFM equation, called the "one-loop" CCFM 
[32,33]. Of course, Q(l — z) is only one of a few ^-dependent functions that have been 
used in the literature [24,25,34]. 

The primary role of the presented here Markovian NLO code would be to serve as a 
powerfull testing device capable of reproducing independently any of the above evolution 
types. In particular, the emulation of the CCFM equation would be possible, if the 
collinear evolution is supplied with the transverse momenta and the "non-Sudakov form 
factor" is added. Another area of application of such a Markovian MC would be in 
performing fits of the PDFs to the F 2 data. We have demonstrated recently [35] that 
using MC codes for the purpose of fitting is feasible with the modern computer power. 
Finally, for the simulation of the final-state parton shower, the Markovian algorithm 
without constraints is directly applicable. 

Let us digress here that the standard evolution kernels of the collinear factorisation, 
either LO or NLO, are used in the normal integrated form, i.e. as functions of the x and 
t variables only. This is of course the well-established and extremely successful DGLAP 
approach. Nonetheless, for more exclusive observables at the LHC, it will be necessary 
to include the effects of transverse momenta beyond the collinear approximation. Some 
attempts in this direction have already been made, for example by the above mentioned 
MC@NLO [11,12], the GR@PPA project [36,37] or the CCFM-based [29] CASCADE 
project [38]. 

Another interesting approach is based on the idea of "exclusive" DGLAP kernels in 
which the internal degrees of freedom are not integrated out analytically, but instead 
simulated (i.e. integrated) by the MC parton shower itself [39]. 

Finally, let us comment on yet another class of MC algorithms solving the evolution 
equations - the veto algorithms [40-42] . They are based on a Markovian evolution with 
additional internal pseudo-emission loop. Thanks to this feature it is enough to compute 
an approximate Sudakov form factor, instead of the exact one, for each of the generated 
events and the algorithm becomes faster. On the negative side, the average multiplicity of 
the generated partons is much higher. We did not use the veto scheme in this work. The 
main reason for this is that, as discussed earlier, we consider Markovian- type algorithms 
primarily as a testing device for the constrained Markovian algorithms. In these latter 
cases the veto algorithm does not apply and the Sudakov form factors must be calculated. 
For that reason we preferred to have the Sudakov form factors implemented (and mutually 
tested) in both approaches. 

This paper is organized as follows. In the next Section we introduce all necessary 
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notation, briefly present the general formalism of the Markovian MC solutions of the 
evolution equations and define two particular evolution schemes that we are interested 
in. In the following two Sections we describe in detail the MC algorithms that solve the 
above two evolutions by means of the collinear Markovian parton shower. We present two 
algorithms, main and auxiliary (mostly for tests), for each of the schemes. As the simpler 
ones we discuss in detail the Q(l — z)-type algorithms, whereas in the description of the 
fey-dependent ones we focus mostly on the modifications with respect to the Q(l — z)- 
type algorithms. In Section 5 we present numerical results: first we discuss the choice 
of the counter term and input parameters, then we briefly describe a variety of technical 
tests that we performed in order to establish the technical precision of the EvolFMC code, 
and finally we compare and discuss various types of evolution. The Summary Section 
concludes the paper. 



2 General framework 

In this paper we will follow the general formulation of the evolution equation and its 
solutions in terms of the Markovian MC presented in Ref. [42]. In order to establish the 
notation let us recall here basic definitions and formulas. For more complete presentation 
we refer the reader to the Ref. [42]. We write the evolution equation in the compact 
matrix notation 

d t B(t) = K(t) D(t) (1) 

and in the explicit representation 



d t D f (t,x) = ^ I dw% 
f Jx 



ff>(t,x,w)D f r(t,w), (2) 



where Df(t,w) denotes the parton density function of the parton / and 3Cfp(t,x,w) 
denotes the generalized evolution kernel. Note that contrary to the DGLAP case, it 
depends on two variables, x and w, instead of their ratio z = x/w only. The kernel is 
then decomposed into real and virtual parts 

%ff(t,x,w) = Xff,(t, x, w) + %ff,(t, x, w), r X V fj,{t ) x 1 w) = —5ff/5 x = w % V ff(t,x). (3) 

The Markovian solution of eq. Q is conveniently expressed with the help of the operator 
E defined as 

ED(t)= / dx y^x D f (t,x), i.e. E f (x)=x. (4) 
Jo f 

This operator in a formal way shows that we use an "unconstrained" evolution - all 
values of x and / will be generated and that the evolution will be done for momentum 
distributions scD(t). The solution of eq. and the master formula for the Markovian 
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MC is then 

ED(t) 



to 



dt\ ( / dt 



dt. 



t 2 



... j dtN\ 'EK R (t N )G- K v(t N , tjv-i) + EG K v(t, tjv_i)<5 tjv=t 

JtN-l L 



(5) 



xK R (t 2 )G K v-(t 2 ,ti) + EG K v(t,ti)(5, 



t 2 =t 



x K R (t i ) G K v (t ! , t ) + EG K v (t , t ) ^ =t ) D (t 



where the diagonal matrix Grv is the solution of the evolution equation with the virtual 
kernel % v ffl {t, x,w) parametrized in terms of the Sudakov form factor <$>f(t,t'\x): 



{G K v(t,t')} fr (x,w) = 5 ff ,5 x=w e -WH $ f (t,t'\x) = f 

Jt> 



dt" W ff (t",x). (6) 



The above formula (J5| follows from the iteration of the momentum conservation principle 

t 

J dti |EK R (tj)G K v(ti, U-i) + EiG K v(t, ti-i)5t i= t | = E. (7) 

ti—i 

Eq. (m) defines also the probabilities of the single step forward in U, fa and Xi variables: 

t 

1 = — / dti (eK^^IGkv^,^!) +EG K v(Mi-i)5 ti=t ) 



X 



J2 dti^fifi-iituU-ilxi-!, j ^ 



X 



(8) 
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where the virtual Sudakov form factor is expressed in terms of the real emission part of 
the evolution kernel 



t, 



U Xi-i 



f J J f 

» U-i Ji 

(9) 

In the case of a weighted algorithm with the global correcting weight, one uses the simpli- 
fied kernel %j.j-._ i (ti,x i ,x i ^i) — > OCf.f. %h a^-i)- This simplification is compensated 
by the global weight 



w (n) = e */„(t,tn[*n)-*/ B (t,tn|*n) | J~J ^A/i-l^' ^ ^ij g*/^-, (t^-ll^-l)-*/,,., 

\i=l ^/i-i^'^i'^-l) 

(10) 

where n denotes the number of emissions generated within a given MC event, and the 
form factor Qf.^ti, is constructed from 3Cj?y. (t i( Xj, in complete analogy 

to eq. ([9]) . The last quantity to be defined here is the exact shape of the kernels, including 
the definition of the argument of the coupling constant in terms of the t and x variables. 
Following Ref. [42] we will discuss two schemes of modified-DGLAP type, denoted in 
Ref. [42] as (B') and (C). The novelty with respect to Ref. [42] is that we will perform 
the calculation and construct the Markovian MC code at the NLO level, whereas in 
Ref. [42] only the LO case has been discussed. To be specific, the schemes are defined as 

x"Xpf \t, x, w) = 



= ( ftMo( ' n( 2 1 ; z)+t) 2zPPf\z) + ( ajV£ ° (ln( 2 1 7r " Z) +t) ) 2 2zPf}p B \z ) y i 
x%f,f \t, x, w) = 



'w—x>Xe~ t 1 
(11) 



where z = x/w and A is the cut-off on the argument of the coupling constant. Note that 
A is not an infinitesimal IR cut-off, as used in the DGLAP case. On the contrary, A is 
finite, typically of the order of a few GeV. Note also that the factor 2 in front of the 
kernels is due to our definition of the evolution time, which in the DGLAP case is chosen 
as t = InQ rather than t = lnQ 2 . 

The NLO parts of the kernels, Pfi^p B (z) and Ppf^ C (z), consist of the universal 
part Pf,}\z) [16,17] and the evolution scheme dependent counter terms AP^j- (z) and 
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APpf )c '(z) 



f 



(12) 



These counter terms are necessary to remove the double counting introduced by the shift 
of the arguments of the coupling constants. There is some freedom in defining these 
counter terms. The algorithms presented here work for any (reasonable) choice of the 
counter terms. For further details on the choice we used in this work we refer the reader 
to Section [5TT1 

We will frequently be using also the representation of the universal parts of the kernels 
based on their structure in the z-variable 
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zP?f\z) = ^-JffAfjiz) + Fpj(z). (13) 

The functions Pp^(z), Aff(z) and Fff(z), written in the notation used here, are collected 
in Ref. [19]. 

The coupling constant at the NLO level has the standard form 

0tL0{t) 



zP flf (z) - ^—d rf A ff + b flf {z 



A^-hiAo)' 

. . / . , /Mn(2t-21nAoA 
a N Lo{t) = a LO (t) I 1 - a LO (t) 1 -J ■ (14) 

We close this general introduction with a brief explanation why we state that the 
argument m(a;.j_i — Xj) + £j can be regarded as kf of the emitted real parton with the 
four momentum k{. It follows immediately from the kinematical mapping of the evolution 
variables into four momenta, provided that the evolution time is identified with the rapid- 
ity of the emitted parton and the x-variable, as usual, with the light-cone plus component 
of the virtual parton 

U = hx(2E h ) - J In %, kf = 2E h (x i _ 1 - Xi ) (15) 
2 fc i 

where Eh is an arbitrary reference energy of the incoming hadron. As a consequence of 
the maslessness of ki we obtain 

kf = Vk+k- = e u (xi-i-Xi). (16) 

Let us now proceed with the description of the novel NLO algorithms for the two 
schemes. 
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3 Markovian algorithm for scheme (B') 



In this section we will present two schemes of solving the evolution B' in a Markovian way 
at the NLO level. We begin with the efficient one and then we present the other scheme, 
devised mostly for the testing purposes. 

3.1 Main algorithm 

The main algorithm is based on the simplified LO DGLAP kernel in which the coupling 
constant is used in the NLO approximation. Specifically we follow the eq. (3.21) of 
Ref. [30], see also [43], and we extend it to the NLO case 

xxfp (t,«,„) ^t + Mi-')) ^,,,^ (17) 

ztffHz) =-L-(S rf Af) + maxF$(3) + M,,,). (18) 

L Z % 

We remind the reader that z = x/w. Note that the singular factor 1/(1 — z) is artificially 
introduced into the F-part in order to achieve analytical integrability of the formula, see 
Ref. [30] for more discussion. The constant My / is defined as 

fo, ifP$°W0, 
I ^ lf P f'f ( z ) = °> 

where rj is a dummy technical parameter. The reason behind introduction of Myf is 
very simple - we want to remove all zeroes in the LO transition matrix, because some of 
them might become non-zero at the NLO level and cause infinite weights. The constant 
i] is therefore added to all kernels that are zero at the LO level. Of course this is purely 
technical, dummy, operation, later on corrected by means of a proper rejection weight. 

The corresponding Sudakov form factor, necessary to generate the time t iy is then 
defined as 

$f (Mi-iki-i) = / l dt ! d(^)^x%ff\t,x,w) 

= - dt dua NLO {t - u){A { ?> + V(maxF}?J(z) + My A) 

* Ju-i Jo z 



rU 
J U-i 



f 

1 dt(p{t,t- t x ) - p(t,0)) (Afj + J2 m ^( F f'f( z ) + M f'f) 

f 



-m) - CiU-MAf} + ^(m f F)^) + My f )), (20) 

/' 
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where 



pit, u) = ^ J du a NLO (t - u) 

If 2 ( I fa\n{2G{t,u)) \ 

nj U p G(t,u) V 2 (3 Q 2 G{t,u) ) 

VPo fa G(t,u)J 



(3 3 G(t,u) 

G(t,u) =t-u-lnAo, u= -ln(l-2;), t A = lnA, (21) 



and 

\2 



C(t) = a 00 (ln(t-lnA )) z + (a 10 t + an) In (t - In A ) + a 20 t, (22) 



° 00 -2/5 3 ' 



ai0 "/V 



an 



1 A 
2 

2 A, 2 In A - Pi In 2 - ft 



A) 



3 



A + 2/j 2 (t A -lnA ) ] 2 / 3 2 (t A -lnA )+ft(l + ln2) 

020 " /5o 3 (t A -lnA ) ln(tA - ln Ao) ^(t A -lnA ) ' 



Let us note that the coefficients need to be calculated only once during initialization 
of the algorithm. The procedure of inverting the function ((ti), necessary to generate the 
time ti, has to be done numerically. 

In the next step we generate the flavor index fa based on the probability 

i/i-i + maXz F fil-i ( z ) + M hh-i \- -t 
">•- WJU^M -4^ /i _ i+E/l (max,F<l iW + M /i/i _ 1 )' Y" ~ ' 

(23) 

where 

(24) 

and 

dt^i^t^x^) = Y^d u $% f ._ l (t i ,t i - X \x i - 1 ). (25) 

Si 
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As the last variable we generate Zj. The normalized density distribution dzip B> \zi) is 
given by 

dzp B '(z) =dz aNL °( ti + ~ Qi-^i^e-'' f f 1 dz aNL °^ ti + ~ Z ^ Ql ~ z>Xe ~ 



7T 1 — Zi \J 7T 1 — Z 

=dv,i aNLO tl — — Ot i -t A>Ui >o (p(ti, U - t\) - p(ti, 0) ) 

=dv,id Ui p(ti, Ui)6 t ._ tA >« i > - £a) - 0) j . (26) 

In the actual generation of we use the method of inverse cumulative. The p(ti, Uj) 
function has to be inverted numerically with respect to the ui variable. 

The last part of the algorithm to be discussed here is the correcting weight, defined 
in eq. (10), compensating the simplification of the kernel done in eq. (17). The most 
complicated part of this weight is related to the exact Sudakov form factor. It has the 
form of the double integral. Numerical evaluation of such a double integral would signif- 
icantly slow down the algorithm. Therefore it is essential to perform at least one of the 
integrations analytically. In the following we will show how this can be done in the NLO 
case. 

Let us define the full Sudakov form factor $j (ti, of the B' evolution 

$f{t i ,t i - 1 \xi- 1 )= f dt I d(-^-}^xXfp(t,x,Xi- X ) 

•J ti 1 *-* 1 j.f 



U-i 



f 



+ + ) h ^-„ (27) 



and let us write down the desired virtual component of the weight 

Af {ti,ti-x\xi-x) =&f(U, ti-i\xi-i) - $f (Mi-iki-i) 

J ti JO f f 

+ { as{t+ ^- z)) )\zlfp B \z)\ (28) 

The integral over dz for general form of the kernel cannot be performed analytically even 
at the LO level (cf. Ref. [30]). However, as in the LO case of Ref. [30], the dt integral 
can be done analytically also for the NLO case. The calculation looks as follows. We 
introduce the usual variable u = — ln(l — z) and then change order of integration. The 
resulting integral can be expressed as a sum of two integrals over the two regions of the 
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u 








ti - t x 










* Jaw 














ii 









Figure 1: The integration domain in the tu space for the case B'. Regions I and II correspond 



to two integrals in eq. (29) 



tu space shown in Fig. 1 L 



pU pi pU pt—tx 

/ dt dz6i^ z> \ e -t — dt due~ u 

Jti-i JO Jti-i Jo 

= / due- 11 I dt+ due- 11 I dt. (29) 

Jti_i—t\ Ju+tx Jo Jti—i 

The dt integral, which depends on the coupling constant only, can be done analytically 
with the help of the integrals 

1 If 

— Ji(t,u) =— / dta NLO (t — u) = —p(t,u) (30) 

71 7C J 

and 



1 1 f 

-J 2 (t,w) = — / dta 2 NLO {t-u) 

(-2 /3 2 (t -u- hiA ) + ft In 2 (t - u - lnA )) : 



7T 

dt 



/? 6 (t-M-lnA ) 



G 3 (t,u] 



\n{2G{t,u)) 



O-20 Q30 a 40 

nil j- \ ~rul \ ' 



G 3 (t,u) G 2 (t,u) G(t,u) 



1 The similar change of the integration order was also exploited by the authors of HERWIG MC [9,44]. 
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where 



IP! 2(31 2 01 Pi 4 



Collecting together the above results we can rewrite eq. (|28|) as 

if 
\f 



| (J 2 (tj,u) - J 2 (u + t X} u)) R{l)B' ( s 
+ 47r 2 ^/'/ W 



^ o f 
+ ^ u) -Jf^ u)) P^ B '(z)) (32) 



where z = 1 — e u . Note the cancellation of the leading singularity, A?J/ (1 — z), in the 
above formula due to 



(33) 



Combining together the real and virtual components we arrive at the final formula for the 
global weight of eq. (10) adopted to the case of the scheme B' 

-R(B') 



w 



fifi- 



n 

i=l 



(34) 



3.2 Auxiliary algorithm 

The main algorithm described in previous section is a fairly complicated one, especially 
due to the presence of numerical integrations and numerical inversions of various functions. 
Therefore, it is obligatory to devise a way of testing it down to the sub-per-mille precision 
level. We have not found any non-Monte-Carlo program that solves the modified-DGLAP- 
type evolutions at the NLO level and therefore we constructed another, independent 
Monte Carlo algorithm for the purpose of cross-checks. The algorithm is less efficient but 
at the same time it is simpler. 

The algorithm is closely based on the LO algorithm described in Ref. [30]. The entire 
NLO correction is introduced as a weight. We will only briefly sketch it here and we refer 
the interested reader to Ref. [30] for details. 
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The algorithm is based on the simplified kernel of the form 

X X?pV(t, X , W ) ^ t + l < 1 -^ 2zPff'\z)e^ z>Xe - t . (35) 



As compared to simplified kernel (17) of the previous algorithm in this one the coupling 



constant is taken in the LO approximation. 



In complete analogy to eq. (|20|) the simplified Sudakov form factor reads 

•£j l) 



Jti-! Jo va^-iV 

= - dt dua LO (t - u)(Af) + YVmaxFf?)^) + M ff )) 

71 Jt i - 1 Jo f z ' 

dt(p LO {t,t- t x ) - p LO (t, 0)) (Af) + ^(maxFp}(z) + M ff ) 



--(Clo(U) - Oo(ti-i))(4/ + ^(maxF^(z) + M ff )), (36) 

/' 



where 



1 f 2 

PLo(t,u) = - / duocLo{t-u) = —— ln(2G(t,u)). (37) 

The function Clo(^) reads 

Cio(ti) = |-(t« - IbAo)(1h(*< - lnA ) - 1) (38) 
Po 

In order to generate t{ the function Clo{U) must be inverted. This inversion can be 
done analytically. However, in the actual Monte Carlo, for the purpose of tests we have 
implemented also the numerical inversion procedure, similar as in the case of the first 
algorithm (there is no visible computing time overhead related to this numerical inversion). 



The generation of the flavor index fa is based on the probability identical to eq. ( 23 ) 



, 2) = ^Mh) = + max ~ 41. M + m (39) 



It is the case because the coupling constants cancel in both eqs. (23) and (39). The 



generation of the z- variable is based on the LO version of eq. (26) 

dz- B ' {2) (z) =dz aL °( tl + ~ Qi-^e- 4 - f f 1 dz aLo((t% + - Z ^ &l ~ z > Xe ~ h \ ~* 

" 7T 1-Z { \J 7T l-Z ) 

=duid Ui pLo(ti,Ui)Bt i -t x >u i >o(pLo(U,U-tx) - Plo(U,0)j ■ (40) 
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Contrary to the previous algorithm, now the function Plo(^ m ) can be inverted analyti- 
cally. However, for the purpose of testing various components of the program, we have 
implemented also an option of the numerical inversion. 

As indicated earlier, the novelty of this algorithm, i.e. the NLO effect, is hidden in the 
global weight. The most complicated part of the weight is the virtual component built 
out of the Sudakov form factors. For the calculation of the exact form factor we can use 
the results derived for the previous algorithm, and the whole virtual part of the weight 
follows from eq. ( 28 ) 



' I 27T / ' / 

/' 



J to JO ft \ 



1/ 

f 



aNLo(t + hi(l - z))\ 2 R(i)B' ( s oi LO (t + ln(l - z)) - R{0) 



+ C 2 ; ' ) ™,>r w - — 2 7 ^ qr'wj. <«) 

with the final result 



rU—t\ / 

Jti-X—t\ fl V 



2tt lZj 



(jlM - Ji(ti + t A ,tQ) fl(o), v {MUjU) - J 2 {u + t x ,u)) R(i)B> 

+ 2tt ^ [Z)+ 4tt 2 ^ 



+ f «2,v( {J ''' >{h -" ] • / ^ / -" )) r;ff , 

^0 V 27T 



2 



| (Jl(Vu) ~ Ji(tj_i,M)) fl(0) (J2(U,U) ~ J 2 (U-1,U)) p R(l)B' ( A 

+ 2vr / ' / 1 J 4vr 2 / ' / 1 J / ' 1 J 

where 2 = 1 — e~" and 

-J?°(t, u)=- ( dta LO {t -u) = -pLo(t, u). (43) 

7T 7T J 

We conclude this section by presenting the complete formula for the global weight in 
this algorithm: 



w {n)(B'){2) = e 



-A^(t,t M "Q I ** ^ c- A ^^-^) . (44) 



4 Markovian algorithm for scheme (C) 

Having completed the presentation of the algorithms for the scheme B' we proceed now 
to the most important scheme C. As discussed in Ref. [42], it can be interpreted as 
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evolution in the rapidity variable with kr as the argument of the coupling constant and it 
is of great physical importance. The NLO algorithms for the scheme C are quite similar 
to the algorithms for the B' scheme. Therefore in the following sections we will skip some 
of the details and concentrate on the differences with respect to the scheme B'. As before 
we begin with the efficient one and then we present the other algorithm used for tests. 



4.1 Main algorithm 

The main algorithm is based on the simplified kernel similar to the one used in the scheme 
B' (i.e. the LO kernel with the NLO coupling constant): 



xXff >(t,x,w) =^i^ _ v -ZzPff \z)0 { i- z)w> xe-^ (45) 

zPf^p (z) =zPfp (z) = (S ff Af) f + max F$ (z) + M ff ). (46) 

Note that now the kernel depends truly on both the x and w variables (through the 9- 
function). In the B' case it depended only on the ratio z = x/w. The simplified Sudakov 
form factor, necessary to generate time ti, is then defined as 



$f {U,U-i\xi-i) = jf dt jf d [j—^j ^xX^'\t,x,x^ x ) 
= — / dt duaNLojt + \nxj-i — u) [AfJ + y](maxF^(z) + M// 

= I dt(p{t + lnxi-i, t + lnx 4 _i - t x ) - pit + lna:i_i, 0)) (a { °} + ^(maxFj?}(«) +M ff ) 

Jti—l ' /' 

= (C{U + hxx^x) - C(ti-i + lnzi-i)) (4/ + £(maxF$(z) + M /v )). (47) 



The functions p and ( are the same as in the algorithm B'. The only difference is in the 
shifted t-argument: t — > t + hiXi-i or > equivalently, shifted integration/generation limits 
ti(i-i) — >■ t*(i-i) + lnxj_i. 

The generation of the flavor index fa is done, identically as in the previous algorithms, 
by means of the probabilities pf t defined in eqs. (23) or (39). 

As the last variable we generate Zi using the integrand of the density function dzip c (zi) 
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given by the formula 

, C' ( \ _ , OLNLojU + baXj-1 + ln(l - Zj)) ©(l-^-^Ae-** 

7T 1 - 

/ f 1 d (^NLpiM +ln(l - z)) ©(l-^^Ae-Q ~' 

\Jo * l ~ z 

_, aNLoiti + lnXi^-Ui) 

—CtUi "ti+lnZi_i-t A >iii>0 

7T 

x (p(ti + Inx^ti + lnxj_i - t\) - p(ti + lnxi__i,0)j 

=dUid Ui p(ti + lriXj_i, Wi)Ot i +lna; i _i-t A >Mi>0 

x ^p(tj + lnxi_i,ti + lnsci-i - t x ) - p(U + lnxi_i,0)j . (48) 

As in the case of tj-variable we observe that the whole difference with respect to the case 
B' is in the shift ti — > t i + lnxj_i. 

Let us define now the full Sudakov form factor $y (ti, 

$f' ) (t J ,t,_ 1 |x l _ 1 )= t dt I d [-^-)Y, xX fP^ 

J ti 1 ^ i 1 r/ 



x,Xi- X ) (49) 



(C r ) 

As in the case B' we are able to perform analytically only one of the integrations in <3>^ . 
To this end we rearrange order of integrations and decompose the integrals as follows 

dt / dz6^ z ) Xi _ 1>Xe -t = dt due~ u 

ti-! JO Jti-! Jo 

(/■ti+lnSj-i-tA rti i-ti-i+hxxi-x—tx rti 

/ due~ u / dt+ due~ u / dt 

Jti-i+\nxi-!-tx Ju—lnxi-i+tx Jo J U-\ 

pti+bxxi-i—tx rti 

+ Otx-xnxi-^ti-, / due~ u / dt (50) 

Jt i - 1 +\nx i - 1 -t x Ju-\nXi-!+tx 
"ti+lriXi-i-tx rti rti-i+lnxi-i—tx fU 

due~ u / dt + / due~ u / dt. 

'ti— i+lnxi_i— tx J u— \nxi-\-\-tx J Jti—i 

(51) 

Note that the decomposition changes when t\ — lnx^-i becomes smaller /greater than 
These two cases are depicted in Fig. |2j With the help of the rearrangement (51) we can 
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Figure 2: The integration domain in the tu space for the case C. Regions I (triangle) and 



II (rectangle) correspond to two integrals in eq. (51), respectively. Left frame : the case of 
t\ — lnXj-i < ti-i. Right frame: the case of t\ — lnxj_i ^ ti-i. 



calculate the virtual component 

Af {UJi-^Xi-i) = $f (Mi-iN-i) -$f (^,^_i|xi_i) 



*i— 1 /"l 

(it / 



o 



(1— z)xi-i>\e 



^ a s (t + lnx t _! + ln(l - z )) \ 2 2z P* wc ' ^ 
2n J f f 



ti+lnxi^x-tx 



tj_i+lnxi_i— 



due- u 2z 



2tt 



J 2 (*i + lnxj-i,tt) - J 2 {u + t\,u) D R(i)c'f \ 



4tt 2 



+ lna;j_i<*i_i 







due~ u 2z 

r 



Jxjtj + lnxj-i,u) - Ji(t,-i +jn^_i 1 u) - fi(0 ) p i?(C)/ v\ 

J 2 (tj + lnXj-i,M) - Mtj-i + lnxj-i,u) u R(i)c' ( s X , _ _„ 

47r 2 ^/'Z Wl' 2-1 e • 



(52) 



Finally, the formula for the global weight, corresponding to eq. (34), in the case C reads 

,(n)(C) = R -A£(t,t„K) JJ f ^gj^' ^l e - A 4 , feA-il*«-i) 



(53) 
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4.2 Auxiliary algorithm 

Similarly to the case B', the second algorithm is based on the version of the simplified 
kernel (45) in which the coupling constant is taken in the LO approximation: 

*X*°™&x,w) = aL ° (t + ' n ^ + h(1 - '» 2zP*P (54) 
As in the case B', in analogy to eq. (47), the Sudakov form factor becomes 

* < j®(t i ,t i - 1 \x i ^ 1 ) = P dt [ d(-) Y^xXff )(2) (t, x, x^) 

Jti-! JO ^ x i-l' j, 

= - dt dua LO (t - u) I Ay} + ( max Ffj (z) + M rf 

= / dt(p LO {t + lnxi_i,t + lnxi_i - t A ) - p LO (^ + lnXj_i,0) 
Af} + J2U^FP f (z)+M f 



\ z 

f 



(Clo{U + ln^_i) - Cxo(*i-i +lnx i _ 1 )) (4/ + S( m f* + M f'f 



(55) 



where 



If 2 

pLoit + hi£i_i, u) = - dua LO {t + lnxj_i - u) = — — ln(2G(t + lnxi_i,w)). (56) 

7T 7 Po 

In order to generate ^ the function C,lo{U + hiXj_i) is then inverted either analytically or 
as a test option also numerically. 

The generation of the flavor index is based on the probability p/ 4 identical to eqs. 
pi and (1391) 



The generation of the z- variable is based on the LO analogue of eq. (48) 
a L o(U + lnxj_! + ln(l - z { )) 9(i_ Si ) !Bi _ 1 >Ae-*< 



dzip c ' (2 \zi) =dz r 



7T 1 - ^ 

1 



1 , OlLoiti + hlXj_i + ln(l - z)) Q(l~ z )x^ 1 >\e'^ 
71 1 — Z 

duid Ui p LO (U + lna;i_i,iXi)e t . +ln!C ._ 1 _i A>tti >o 

(plo(U + \nx i _ 1 ,t i + lnxj_i -t x ) - Pw^i + lnx i _i,0)) . (57) 



As in the case B', the whole NLO effect is introduced through the global weight. The 
virtual part of this weight, constructed on the base of the exact Sudakov form factor 
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derived for the previous algorithm, follows from eq. (52) 

A? (2) (ti,ti_i|x^i) = $f {t^U^x^x) - $f (2) (Mi-ik*-i) 

At f 1 a a sr^faNLo{t + lnXi_i + ln(l - z)) m 

= Jt Jo ( 1 - 2 ) x - 1>Ae " t 2^ 2zP/ f '{z) 



f 

^NLo(t + lngj-i + ln(l - z )) \ 2 2 zP R{1)C ' (z) 
2% J * * 



_ c LO (t + lnx t _ 1 + ln(l-,)) 2 ^ (0) 

27T 7 7 

leading to 



Jti-1+lD.Xi-!-^ f, V 27T 

(Ji(ti + lnxj_i,tt) - Ji(M + t A ,w)) b( )/ x , (MU + *Q ~ + t A , tp) k( 

+ 2 7T ^ [Z)+ 4tt 2 ^ 

j-ti-i+Ina^-i— *a 

+ fl^-ina*-!^-! / due~ u 2z 

JO 

E( {Ji°{ti + \nxi^u) - Ji {U_i + lnXj_i,M)) pfl(C) r x 
V 2vr W 

| {Ji{U + lnxj_i,^) - Ji(tj_i + lnxj_i,^)) fl(p), 
+ 2vr 

| (J 2 (tj + \nxj_ u u) - J 2 (Vi + lnxj_i,tt)) p fl(i)cy A _ 1 _ u rc . Q x 
+ ^ 2 ^ /;/ (2H , 2 - 1 e . <ayj 



As a result, the complete formula for the global weight relevant for this algorithm reads 
,(■0(0(2, = e -A^(t M "Q ( X ^ i,Xi,Xi - l] e -±£>*-i\«-S\ . (60 ) 



-R(C) 

fifi — 1 



5 Numerical results 

In this Section we present numerical results obtained with the MC program EvolFMC 
version 2. The first version of this program has been presented in the Ref. [18] for the 
case of standard DGLAP LO evolution. Subsequently, the NLO evolution has been added 
in Ref. [19] and the LO modified-DGLAP evolution in Ref. [30]. The presented here 
EvolFMC v.2 includes three of the described above four algorithms for NLO modified- 
DGLAP evolution (the auxiliary algorithm for the B' evolution is not implemented at the 
moment). In order to accommodate these new evolution schemes, the overall structure 
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of the program has been modified, see Ref. [39] for details. As a result, it was important 
to perform a number of technical comparisons of the new code in order to establish its 
technical precisions. We very briefly describe these tests in the following subsection. 
Later on we present numerical results regarding the actual new evolution schemes. Before 
showing the results we discuss the choice of the counter terms and we list the input 
parameters. 



5.1 Removal of the double counting 

As indicated in Section [2j modification of the argument of the coupling constant in the 
evolution kernels is equivalent to adding some higher order (i.e. NLO and higher) terms. 
Therefore, one has to make sure that there is no double counting with the NLO part of 
the kernel. In the implementation presented here we follow the approach as discussed in 
Ref. [34]. Namely, we take the expansion of the aNLo{t) in the form 

a NLO (t + ln0) = a NLO (t) - ^a 2 NLO (t) ln0 + 0(l/t 3 ), (61) 

where <fi represents any arbitrary change in the argument. Then the extra term in the 
kernels is of the form 

-^{^^)\zPff (z) * -/3 o ln0( a ^°^ +ln ^ ) 2 2<f(,) + 0(l/t% 

(62) 



and consequently the counter terms from eqs. (12) could be defined as follows 



AP*f )B \z) = p ln(l - z)Ppf\z), for scheme B', (63) 
APff C \z) =/? (lnw + ln(l - z))P*f\z), z = x/w, for scheme C\ (64) 



However, a more detailed inspection of eq. (64) reveals that this formula over- subtracts 



the double counting. Namely, in the DGLAP evolution the kernels are functions of "local" 



z- variables only. In eq. (64) there is a corresponding ln(l — z) term. The other term, the 
lnu>, interconnects two emissions, and as such it is of genuinely beyond-DGLAP origin. 
It means that it is absent in the DGLAP kernel and there is no reason to remove it. As 
a result, the counter term in the C case becomes identical to the counter term in the B' 
case. In the following we will present results for both variants in the C case, with the 
understanding that the preferred choice is the 

APff B \z) = AP B f )C '(z) = (3 ln(l - z)P*f\z), (65) 

for both schemes B' and C'. 
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5.2 Input parameters 

The set-up of the EvolFMC code is the same for all the presented results. We use the gluon 
(G) and quark singlet (Q) PDFs with three massless quarks 

d q = J2{ d ^ + d ^)- ( 66 ) 

i 

As the initial distributions of the evolution we take 

D%(x) = 1.908 • x~ 12 {l - x) 5 - , 
D° sea (x) = 0.6733 -x-^il-x) 7 - , 

D° u (x) = 2.187 -x-^il-x) 3 - , 1 j 

DlJx) = 1.230 .x-°- 5 (l-xr 

and 

D° u (x) = Dl ai (x)+ 1 -Dl a (x), 

D° d (x) = DlJx) + 1 -Dl a (x), (6g) 

D° s (x) = Dl{x) = D%x) = D°(x) = \Dl a {xl 
Dl(x) = DUx) + D° Uv Jx) + DlJx). 

The QCD constant A = 0.2457, the cut-off A = 1 GeV, Nf = 3 and the dummy parameter 

T) = 0.1. 



5.3 Technical tests 

We have performed three different sets of the technical comparisons of the code EvolFMC 
v.2: 

1. With the semianalytical code APCheb40 [42,45] based on the expansion in Chebyshev 
polynomials. 

2. With the previous version of the EvolFMC code: the EvolFMC v. 1, which was ex- 
tensively tested in the past. 

3. Between different algorithms within the EvolFMC v.2. It is the most important test 
of the new NLO modified-DGLAP evolutions, which are not available in any other 
code. 

The overall conclusion of all the tests is that the technical precision of the program 
EvolFMC v. 2 is at least 5 x 10~ 4 (half of a per mille). For the details of the tests we refer 
the reader to Ref. [39]. 
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5.4 Comparison of evolutions 

Having established the technical precision of the EvolFMC code let us now proceed with 
the comparison of the various evolutions discussed earlier. Before getting into details let 
us remind the reader that in this paper for all of the results (i.e. for all of the evolutions) 
we use the same parameter setup: the same initial distributions at Q = 1 GeV and the 
same A . In a more realistic study one should perform fits to the experimental data 
for each of the evolutions separately and then use different initial distributions for each 
evolution. 

We organize the numerical comparisons in the following way: We begin by showing 
the "reference" result, i.e. the standard DGLAP in the LO and NLO approximations 
(Fig. [3]). Next, we show the new results for the B'-type and C'-type evolutions. These 
two plots (Figs. [4] and [5]) are the main numerical results of this paper, showing the NLO 
corrections in the modified DGLAP evolutions. The next two plots (Figs. [6] and [7]) are 
of technical character and show in a more detail certain aspects of the C'-type evolution. 
We conclude the section by comparing in a common plot all three evolutions in the NLO 
approximation (Fig. [8]). 

5.4.1 NLO versus LO 

The "reference" DGLAP evolution. In the Fig. [3] we show the case of the standard 
DGLAP. We show these well known results because DGLAP will serve us as a reference 
point in discussing the modified evolutions of the B'- and C'-type. We present the gluon 
and quark momentum distributions in LO and NLO approximations as well as their ratios. 
Three evolution time limits are shown: 10, 100 and 1000 GeV. The characteristic feature 
of the plots is that DGLAP NLO corrections are systematically bigger in the large x 
region and they show the tendency of diverging in the x — > 1 limit. In the small x region 
the NLO corrections are small. Results for the other evolutions will be presented in a 
similar way. 

The B'-type evolution is shown in the Fig. |4j It is one of the two main new numerical 
results of this paper. As compared to the DGLAP case we can notice that the NLO 
corrections are a bit bigger in the small x region, of the size up to 20%. In the large 
x region, on the contrary, the corrections are much smaller and showing less divergent 
behavior. This is in agreement with the general principle of discussed here modifications 
of the DGLAP equation. These modifications are supposed to improve the description of 
the emission of soft partons [23,26,28] and it is the x — > 1 limit which corresponds to the 
limit of only soft emissions. 

The C'-type evolution. In the Fig. [6] we compare the LO and NLO evolutions in the 
case of modified-DGLAP C'-type. We consider this figure as the most important new 
numerical result of this paper. We show the gluon and quark momentum distributions for 
two evolution time limits: 10, 100 and 1000 GeV. We see that at 100 and 1000 GeV the 
NLO corrections in the small-x region are somewhat bigger than in the case B', reaching 
even 30%. In the large x region the NLO corrections seem to be even milder and slightly 
less divergent than in the B' case, showing the same improvement over DGLAP. For the 
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Figure 3: Left frames : the DGLAP-type evolutions in the LO and NLO approximations. Upper 
curves (LO: magenta and NLO: blue): the gluon xDg(x) distr.; lower curves (LO: black and 
NLO: red): the quark xDq(x) distr. Right frames: the ratio of the NLO to LO distributions 
for the gluon (magenta) and quark (black) distributions. Top frames: the evolution up to 10 
GeV. Bottom frames : the evolution up to 100 GeV. 



shorter time of 10 GeV the effects are much smaller. In fact, in this evolution type, due 
to the cut-off hp > A, much less of the evolution happens before 10 GeV, and both the 
LO and NLO curves are close to the initial condition as well as to each other. Note that 
in the case of both DGLAP and B', the evolution at 10 GeV is already well developed. In 
order to reduce the evolution and get closer to the initial condition, in the DGLAP and 
B' cases the evolution time must be much shorter, below 2 GeV at least. 
Technical details related to C'-type evolution: 

The C'-type evolution with the disfavoured counter term AP from eq. (64)- For the sake of 
comparison, in the Fig.[6j we present also the other, disfavored, choice of the counter term 
in the C evolution, given in eq. (64). The plots clearly confirm that the NLO corrections 
are much bigger, and, in addition, strongly divergent in the small x limit. 
The C'-type evolution with no counter term. As the last exercise, in Fig. [7| we completely 
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Figure 4: Left frames : the modified-DGLAP B'-type evolutions in the LO and NLO approxi- 
mations. Upper curves (LO: magenta and NLO: blue): the gluon xDg(x) distr. ; lower curves 
(LO: black and NLO: red): the quark xDq(x) distr. Right frames: the ratio of the NLO to 
LO distributions for the gluon (magenta) and quark (black) distributions. Top frames: the 
evolution up to 10 GeV. Bottom frames : the evolution up to 100 GeV. 

switched off the counter term in the C'-type evolution. This way we try to understand 
better what actually causes the changes in the shape of the NLO corrections in the 
modified schemes: is it the genuine change of the argument of the coupling constant or 
rather the cut-off A? In a crude approximation the counter term can be regarded as an 
estimate of the size of the pure effect of shifting the argument in the coupling constant. 
As one can see, the shapes in the right hand side plots in Fig. [7] have changed significantly 
in the large x region, and have remained similar in the small x region, as compared to 
the complete C plots from Fig. [5] This demonstrates that indeed it is the change of the 
argument that drives the effect in the region of soft emission. On the other hand, from 
the comparison of Figs. [5] and [7] to the DGLAP evolution, Fig. |3j we inferr that the higher 
order terms, resumed in the coupling constant, contribute as much as the counter term. 
Remark on the small x limit of the C'-evolution. As we have mentioned, C evolution 
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Figure 5: Left frames : the modified-DGLAP C'-type evolutions in the LO and NLO approxi- 
mations. Upper curves (LO: magenta and NLO: blue): the gluon xDg(x) distr.; lower curves 
(LO: black and NLO: red): the quark xDq(x) distr. Right frames: the ratio of the NLO to 
LO distributions for the gluon (magenta) and quark (black) distributions. Top frames: the 
evolution up to 10 GeV. Central frames : the evolution up to 100 GeV. Bottom frames : the 
evolution up to 1000 GeV. 

is motivated by the CCFM equation. This equation treats the x — > limit in a better 
way than the DGLAP equation, in a sense interpolating between DGLAP and BFKL 
equations. However, apart from the modifications in the coupling constant and the angular 
ordering, which we have incorporated in the C scheme, there is one more ingredient 
missing in C - the 'non-Sudakov' form factor. This non-local form factor strongly modifies 
the small x evolution. It depends on the transverse momenta of all the emitted partons. 
It is in principle not difficult to include such an effect into the Monte Carlo evolution of 
the C'-type by means of the rejection mechanism, provided the transverse momenta are 
properly generated within the cascade, and we plan to do it in the future [39]. 

Let us summarize the comparisons: 1. The modifications of the argument of the 
coupling constant decrease the size of the NLO corrections and reduce the divergences 
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Figure 6: Left frames : the modified-DGLAP C'-type evolutions in the LO and NLO approx- 
imations with the modified disfavoured (C'-type) counter term. Upper curves (LO: magenta 
and NLO: blue): the gluon xDq(x) distr.; lower curves (LO: black and NLO: red): the quark 
xDq(x) distr. Right frames: the ratio of the NLO to LO distributions for the gluon (magenta) 
and quark (black) distributions. Top frames: the evolution up to 10 GeV. Central frames : the 
evolution up to 100 GeV. Bottom frames : the evolution up to 1000 GeV. 

in the large x region. 2. The NLO corrections in schemes B' and C are significant, up 
to 30% in small x region and must be included in the MC parton showers for the LHC. 
3. One must remember, however, that these comparisons are to some degree artificial 
because, as discussed earlier, we use the same input PDF distributions for all evolution 
types. In a more complete study each of the evolutions should be separately fitted to 
the experimental data and obtained this way different input PDF distributions should be 
used. 
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Figure 7: Left frames : the modified-DGLAP C'-type evolutions in the LO and NLO approx- 
imations without the counter term. Upper curves (LO: magenta and NLO: blue): the gluon 
xDq(x) distr. ; lower curves (LO: black and NLO: red): the quark xDq(x) distr. Right frames: 
the ratio of the NLO to LO distributions for the gluon (magenta) and quark (black) distribu- 
tions. Top frames: the evolution up to 10 GeV. Central frames : the evolution up to 100 GeV. 
Bottom frames : the evolution up to 1000 GeV. 

5.4.2 Different evolutions 

As the last exercise we compare the three analyzed types of the evolution. In the Fig. 
[8] we show simultaneously all three evolutions in the NLO approximation. It is clearly 
visible that the difference between DGLAP and B'-type evolution (with a(Q(l — z))) is 
rather small and of the quantitative form. On the contrary, the C'-type evolution (with 
a{hr)) looks very different, both in the shape and in the magnitude. There is a visible 
flattening of the C'-type distributions around the value of x ~ X/Q caused by the cut-off 
A. This follows from the condition Xj_i(l — Zj) > X/Q which stops any cascade as soon 
as the Xi variable falls below the X/Q threshold. In particular no evolution at all can 
develop for any x < X/Q. 



26 




Iog 10 (x) log l0 (x) 



Figure 8: Comparison of three NLO evolutions: DGLAP (magenta), B'-type (blue) and C- 
type (black). Left frames : the gluon xDg(x) distr. Right frames: the quark xDq(x) distr. 
Top frames: the evolution up to 10 GeV. Central frames : the evolution up to 100 GeV. 
Bottom frames : the evolution up to 1000 GeV. 

6 Summary and outlook 

In this paper we have presented a series of Markovian MC algorithms that solve the 
QCD evolutions of the modified-DGLAP type in the NLO approximation. One of the 
two discussed modifications of the DGLAP evolution is of high practical importance. In 
this evolution, as the argument of the coupling constant the transverse momentum of the 
emitted parton is used, and the evolution time is identified with the rapidity variable. 
Such a modification is known to describe better the emission of soft partons. It can serve 
as a first step towards incorporating the complete CCFM effects into the evolution as well. 
In this paper we have called this scenario the C'-type evolution. The other scenario, called 
throughout the paper the B'-type evolution, uses Q(l — z) as the argument of the coupling 
constant, modelling the "one- loop" CCFM evolution equation. It is one of many possible 
z-dependent modifications of the argument discussed in the literature. Proper counter 
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terms have been added to the evolution kernels in order to remove double counting at the 
NLO level. 

The algorithms have been implemented into the new version of the MC program 
EvolFMC and extensively tested by comparisons with the semianalytical code APCheb40, 
with the previous version of EvolFMC and by comparisons of independent algorithms 
within the new version of EvolFMC itself. As the overall conclusion of the tests we claim 
the technical precision of the EvolFMC to be at least 5 x 10 -4 . 

The comparison of the modified-DGLAP evolutions at both LO and NLO level shows 
that the NLO corrections are in general smaller in the modified evolutions and more 
importantly, the divergent behavior of the NLO corrections in the large x region is limited, 
as expected. The only exception is the very low x region where the NLO corrections in 
the modified schemes are larger. This is however the region where the DGLAP equation 
becomes less accurate anyway. 

Quantitatively, the NLO corrections are relatively modest: in the B' case they are 
small, of the order of up to 20% of the LO terms. For the ^-dependent evolution (C) 
they are, in most of the parameter space, of the order of 10%, but in some limited 
regions of small x they grow up to 30%. These results, on the one hand, show that 
the NLO contributions are numerically significant and should be taken into account in 
the construction of the parton shower MC for the LHC experiments. This is especially 
important in the case of the physically well motivated /cy-dependent evolution. On the 
other hand, the convergence of the QCD perturbative expansion looks reasonably well, 
better than in the DGLAP case. 

The main limitations of the presented in this paper MC algorithms are: missing effects 
due to the non-zero masses of the quarks in all of the discussed types of evolution and 
lack of the dedicated fits to the data for the modified-DGLAP type evolutions. Another 
interesting development line of the modified-DGLAP type evolutions would be to include 
a non-perturbative parametrization of the behavior of the coupling constant below the 
Landau pole. We hope to address some of these issues in the future. 
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